

clear all
cap log close
set more off
set matsize 1600
set mem 100m

global ROOT ""
global data "$ROOT/Data"
global dofiles "$ROOT/DoFiles"
global tables "$ROOT/Tables"
global figures "$ROOT/Figures"

u "$risk_data\analysisfinal.dta", clear

*-------------------------------------------------------------------
*Deciding on the zones and villages
*-------------------------------------------------------------------
global ourcluster "zone_o"

*Village variable for the village fixed effects
*vill1 is for the village in the first wave
local village "vill1"
 
/* First, generating interaction dummies between the shocks and women's group treatment variable */
qui {

gen wgonly_o = 1 if treatment_o==2
replace wgonly_o = 0 if treatment_o==1 | treatment_o==3 | treatment_o==4


*Shock variables
foreach y of varlist dcrop dintensity {
gen `y'_treat1 = `y'*wgonly_o

}

/* next creating interaction terms between the shock variable and an indicator for the control group*/
gen control = wg_o ==0
replace control = . if wg_o==.

foreach y of varlist dcrop dintensity {
gen `y'_control = `y'*control

}

}
*closing quietly above


*-------------------------------------------------------
* Removing outlying values for savings, livestock
*-------------------------------------------------------
replace large_an = . if large_an>100
replace small_an = . if small_an>120
replace savings = . if savings>500000

xtset num wave
foreach y of varlist large_an small_an savings {
replace d`y' = d.`y'
}

*--------------------------------------------------------
*Building some variables that were not created before
*--------------------------------------------------------
gen d12_18=dfem12_18+dmale12_18
gen dmt18=dfemmt18+dmalemt18

cap program drop banyy
program define banyy, eclass
args v1
sum `v1' if wg_o==0 
local mean = r(mean)
ereturn scalar mean = `mean'
end


/* Next for the regression. The regression we run is as follows:
dsavings = beta*dshock(ivt) + beta2*dshock(ivt)*wg_o + vt + uivt
It comes from the model of perfect risk sharing with CRRA utility. Test is similar to that used by Mace (1991) and Cochrane (1991)
 */


 

********************************************
 * Specification with village-time dummies
*********************************************
cd "$figures"

global regressors1 " ddays_oct dlt6 dsix_12 d12_18 dmt18"
global outregopts "addstat("Intracluster Correlation", r(rho)) bracket"


// Creating aggregates using the Anderson (2008) method.
*------------------------------
* Creating the summary index
*------------------------------
* Step 1: Standardising the variables to create effect sizes
foreach y of varlist dsavings dlarge_an dsmall_an {
sum `y'
local m = r(mean)
sum `y' if wg_o==0
local s = r(sd)
gen std_`y' = (`y' - `m')/`s'
}

regress std_dsavings std_dlarge_an std_dsmall_an
keep if e(sample)


*Step 2: Creating the weights matrix
cap mata clear
mata
st_view(X=., ., ("std_dsavings", "std_dlarge_an", "std_dsmall_an"))
Y = invsym(X'X)			/* Computing the inverse of the covariance matrix */
b = rowsum(Y)			/* Getting the sum across rows */
c= colsum(b)
/* Weights matrix is now */
weights = b/c
st_matrix("weights", weights)

end

foreach y of numlist 1/3 {
local we`y' = weights[`y',1]
}

* Creating the chats index
gen savings_index= `we1'*std_dsavings+`we2'*std_dlarge_an+`we3'*std_dsmall_an

mat drop weights 

matrix coeff1 = J(3,4,.)
matrix colnames coeff1 = "Savings index" "Savings (10000s MK)" "Large Animals" "Small Animals"
matrix CI1 = J(2,4,.)
matrix colnames CI1 = "Savings index" "Savings (10000s MK) " "Large Animals" "Small Animals"
matrix rownames CI1 = ll95_shock ul95_shock 
matrix CI2=J(2,4,.)
matrix colnames CI2 = "Savings index" "Savings (10000s MK)" "Large Animals" "Small Animals"
matrix rownames CI2 = ll95_shock ul95_shock
matrix CI3=J(2,4,.)
matrix colnames CI3 = "Savings index" "Savings (10000s MK)" "Large Animals" "Small Animals"
matrix rownames CI3 = ll95_shock ul95_shock

//For intensity now
matrix coeff2 = J(3,4,.)
matrix colnames coeff2 = "Savings index" "Savings (10000s MK)" "Large Animals" "Small Animals"
matrix CI12 = J(2,4,.)
matrix colnames CI12 = "Savings index" "Savings (10000s MK)" "Large Animals" "Small Animals"
matrix rownames CI12 = ll95_shock ul95_shock 
matrix CI22=J(2,4,.)
matrix colnames CI22 = "Savings index" "Savings (10000s MK)" "Large Animals" "Small Animals"
matrix rownames CI22 = ll95_shock ul95_shock
matrix CI32=J(2,4,.)
matrix colnames CI32 = "Savings index" "Savings (10000s MK)" "Large Animals" "Small Animals"
matrix rownames CI32 = ll95_shock ul95_shock

//Now for the regression
areg savings_index dcrop dcrop_treat1         $regressors1, absorb(vill1) cluster ($ourcluster)
matrix coeff1[1,1]=e(b)[1,1]
matrix coeff1[2,1]=e(b)[1,2]
boottest dcrop , boot(wild) seed(10101) bootcl($ourcluster)
matrix CI1[1,1]=r(CI)[1,1]
matrix CI1[2,1]=r(CI)[1,2]
boottest dcrop_treat1 , boot(wild) seed(10101) bootcl($ourcluster)
matrix CI2[1,1]=r(CI)[1,1]
matrix CI2[2,1]=r(CI)[1,2]

qui areg savings_index dcrop_control dcrop_treat1  $regressors1, absorb(vill1) cluster ($ourcluster)
matrix coeff1[3,1]=e(b)[1,2]
boottest dcrop_treat1 , boot(wild) seed(10101) bootcl($ourcluster)
matrix CI3[1,1]=r(CI)[1,1]
matrix CI3[2,1]=r(CI)[1,2]

//Incidence...
areg savings_index dintensity dintensity_treat1         $regressors1, absorb(vill1) cluster ($ourcluster)
matrix coeff2[1,1]=e(b)[1,1]
matrix coeff2[2,1]=e(b)[1,2]
boottest dintensity , boot(wild) seed(10101) bootcl($ourcluster)
matrix CI12[1,1]=r(CI)[1,1]
matrix CI12[2,1]=r(CI)[1,2]
boottest dintensity_treat1 , boot(wild) seed(10101) bootcl($ourcluster)
matrix CI22[1,1]=r(CI)[1,1]
matrix CI22[2,1]=r(CI)[1,2]

qui areg savings_index dintensity_control dintensity_treat1  $regressors1, absorb(vill1) cluster ($ourcluster)
matrix coeff2[3,1]=e(b)[1,2]
boottest dintensity_treat1 , boot(wild) seed(10101) bootcl($ourcluster)
matrix CI32[1,1]=r(CI)[1,1]
matrix CI32[2,1]=r(CI)[1,2]

local i=1

gen dtsavings=dsavings/10000

foreach y in tsavings large_an small_an {

local ++ i

**Incidence
areg d`y' dcrop dcrop_treat1         $regressors1, absorb(vill1) cluster ($ourcluster)
matrix coeff1[1,`i']=e(b)[1,1]
matrix coeff1[2,`i']=e(b)[1,2]
boottest dcrop , boot(wild) seed(10101) bootcl($ourcluster)
matrix CI1[1,`i']=r(CI)[1,1]
matrix CI1[2,`i']=r(CI)[1,2]
boottest dcrop_treat1 , boot(wild) seed(10101) bootcl($ourcluster)
matrix CI2[1,`i']=r(CI)[1,1]
matrix CI2[2,`i']=r(CI)[1,2]

qui areg d`y' dcrop_control dcrop_treat1  $regressors1, absorb(vill1) cluster ($ourcluster)
matrix coeff1[3,`i']=e(b)[1,2]
boottest dcrop_treat1 , boot(wild) seed(10101) bootcl($ourcluster)
matrix CI3[1,`i']=r(CI)[1,1]
matrix CI3[2,`i']=r(CI)[1,2]

//Intensity...
areg d`y' dintensity dintensity_treat1         $regressors1, absorb(vill1) cluster ($ourcluster)
matrix coeff2[1,`i']=e(b)[1,1]
matrix coeff2[2,`i']=e(b)[1,2]
boottest dintensity , boot(wild) seed(10101) bootcl($ourcluster)
matrix CI12[1,`i']=r(CI)[1,1]
matrix CI12[2,`i']=r(CI)[1,2]
boottest dintensity_treat1 , boot(wild) seed(10101) bootcl($ourcluster)
matrix CI22[1,`i']=r(CI)[1,1]
matrix CI22[2,`i']=r(CI)[1,2]

qui areg d`y' dintensity_control dintensity_treat1  $regressors1, absorb(vill1) cluster ($ourcluster)
matrix coeff2[3,`i']=e(b)[1,2]
boottest dintensity_treat1 , boot(wild) seed(10101) bootcl($ourcluster)
matrix CI32[1,`i']=r(CI)[1,1]
matrix CI32[2,`i']=r(CI)[1,2]

}

mat ren coeff1 Incidence
mat ren coeff2 Intensity

set scheme s1mono
 
coefplot (matrix(Incidence[1,]), ci(CI1) label(Control)) (matrix(Incidence[3,]), ci(CI3) label(Treatment)) (matrix(Incidence[2,]), ci(CI2) label("Treatment Effect")) ,  xline(0) recast(bar) ciopts(recast(rcap)) citop barwidth(0.25)
graph save "$figures/savings_incidence", replace
graph export "$figures/savings_incidence.pdf", replace

coefplot (matrix(Intensity[1,]), ci(CI12) label(Control)) (matrix(Intensity[3,]), ci(CI32) label(Treatment)) (matrix(Intensity[2,]), ci(CI22) label("Treatment Effect")) ,  xline(0) recast(bar) ciopts(recast(rcap)) citop barwidth(0.25)
graph save "$figures/savings_intensity", replace
graph export "$figures/savings_intensity.pdf", replace


coefplot (matrix(Incidence[1,]), ci(CI1) label(Control)) (matrix(Incidence[3,]), ci(CI3) label(Treatment)) (matrix(Incidence[2,]), ci(CI2) label("Treatment Effect")) || (matrix(Intensity[1,]), ci(CI12) label(Control)) (matrix(Intensity[3,]), ci(CI32) label(Treatment)) (matrix(Intensity[2,]), ci(CI22) label("Treatment Effect")),  xline(0)  ciopts(recast(rcap)) citop graphregion(fcolor(white))
graph save "$figures/savings_combined", replace
graph export "$figures/savings_combined.pdf", replace

cap log close
exit
